function jts = calcAhTInv(invP,etaModel,w)

k = 1/invP.disT(1,2);


T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - etaModel)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - etaModel)));
M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

A = invP.mesh.CURL'*(S\T)*invP.mesh.CURL + M;

hts = cell(invP.nt,1);
ets = cell(invP.nt,1);
jts = cell(invP.nt,1);

rhs = w{end};

hts{end} = A\rhs;
ets{end} = (S\T)*invP.mesh.CURL*hts{end};
jts{end} = -T\ets{end};

for ii = invP.nt-1:-1:1
   rhs = w{ii} + M*hts{ii+1} + invP.mesh.CURL'*((S\T)*St - Tt)*jts{ii+1};
   hts{ii} = A\rhs;
   ets{ii} = (S\T)*invP.mesh.CURL*hts{ii} + (Tt - (S\T)*St)*jts{ii+1};
   jts{ii} = -T\ets{ii} + (T\Tt)*jts{ii+1};
end